Chapter 7 Beta diversity

load("data/data.Rdata")
load("data/beta.Rdata")

7.1 Captive-bred vs wild

load("data/beta_filtered_capwild.Rdata")
load("data/data_host_filtered.Rdata")
wild_pre_samples <- sample_metadata %>% 
  filter(diet!="Post_grass") %>%
  dplyr::select(sample) %>% pull()

genome_counts<- genome_counts %>% 
  column_to_rownames("genome") %>% 
  select(all_of(wild_pre_samples))%>% 
  rownames_to_column("genome")

sample_metadata <- sample_metadata %>% 
  filter(diet!="Post_grass")
beta_q0n <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1)

genome_counts <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0))%>%
  rownames_to_column(., "genome")

genome_tree <- keep.tip(genome_tree, tip=genome_counts$genome)

beta_q1p <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, tree = genome_tree)

genome_counts_filt <- genome_counts[genome_counts$genome %in% rownames(genome_gifts),]
genome_counts_filt <- genome_counts_filt %>%
  remove_rownames() %>% 
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0))%>%
  rownames_to_column(., "genome")

genome_gifts1 <- genome_gifts[rownames(genome_gifts) %in% genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]

dist <- genome_gifts1 %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

beta_q1f <- genome_counts_filt %>%
#  remove_rownames() %>% 
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, dist = dist)

7.1.1 Permanova analysis

set.seed(1234)

Richness

sample_metadata_row<- column_to_rownames(sample_metadata, "sample") 
sample_metadata_row <- sample_metadata_row[labels(beta_q0n$S), ]

betadisper(beta_q0n$S, sample_metadata_row$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.004224 0.0042240 7.8055    999  0.003 **
Residuals 22 0.011905 0.0005412                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Pre_grass  Wild
Pre_grass           0.008
Wild       0.010583      
adonis2(beta_q0n$S ~ diet, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
tinytable_y2dqa26z39giqc2kfisq
term df SumOfSqs R2 statistic p.value
diet 1 0.01526961 0.2489907 7.29391 0.001
Residual 22 0.04605641 0.7510093 NA NA
Total 23 0.06132602 1.0000000 NA NA

Neutral

betadisper(beta_q1n$S, sample_metadata_row$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.006925 0.0069254 0.5404    999  0.491
Residuals 22 0.281957 0.0128162                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Pre_grass  Wild
Pre_grass           0.476
Wild        0.47005      
adonis2(beta_q1n$S ~ diet, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
tinytable_mlxrdgkl0941ux9o1nm6
term df SumOfSqs R2 statistic p.value
diet 1 1.538238 0.2539205 7.487472 0.001
Residual 22 4.519715 0.7460795 NA NA
Total 23 6.057953 1.0000000 NA NA

Phylogenetic

betadisper(beta_q1p$S, sample_metadata_row$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.007671 0.0076708 1.8217    999  0.179
Residuals 22 0.092637 0.0042108                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Pre_grass  Wild
Pre_grass           0.187
Wild        0.19084      
adonis2(beta_q1p$S ~ diet, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
tinytable_3f450bngrbc1njb7kdqe
term df SumOfSqs R2 statistic p.value
diet 1 0.3720264 0.3082818 9.804858 0.001
Residual 22 0.8347476 0.6917182 NA NA
Total 23 1.2067741 1.0000000 NA NA

Functional

betadisper(beta_q1f$S, sample_metadata_row$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.003070 0.0030698 1.7361    999  0.171
Residuals 22 0.038901 0.0017682                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Pre_grass  Wild
Pre_grass           0.176
Wild        0.20119      
adonis2(beta_q1f$S ~ diet,
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
tinytable_glmf4taaiuawxafp2s1w
term df SumOfSqs R2 statistic p.value
diet 1 0.01079711 0.187686 5.083124 0.103
Residual 22 0.04673041 0.812314 NA NA
Total 23 0.05752752 1.000000 NA NA

7.1.2 Beta diversity plots

7.1.2.1 Richness

beta_q0n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(diet) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
    geom_point(size = 4) +
  scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Post_grass" = "Captive-bred grass diet", "Wild" = "Wild"))+
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 12, face = "bold"),
      axis.text = element_text(face = "bold", size = 12),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 12),
      legend.title = element_text(size = 14),
      legend.position = "right", legend.box = "vertical"
    ) +labs(color='Origin')

7.1.2.2 Neutral

7.1.2.3 Phylogenetic

beta_q1p$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(diet) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
    geom_point(size = 4) +
    scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Wild" = "Wild"))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(color='Origin')

7.1.2.4 Functional

beta_q1f$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(diet) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
  geom_point(size = 4) +
  scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Wild" = "Wild"))+
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme_classic() +
  theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(color='Origin')

7.2 Before and after grass

load("data/beta_filtered_pre_post.Rdata")
genome_counts <- genome_counts_filt
post_pre_samples <- sample_metadata %>% 
  filter(diet!="Wild") %>%
  dplyr::select(sample) %>% pull()
genome_counts<- genome_counts %>% 
  column_to_rownames("genome") %>% 
  select(all_of(post_pre_samples))%>% 
  rownames_to_column("genome")
sample_metadata <- sample_metadata %>% 
  filter(diet!="Wild")
beta_q0n <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1)

genome_counts <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0))%>%
  rownames_to_column(., "genome")

genome_tree <- keep.tip(genome_tree, tip=genome_counts$genome)

beta_q1p <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, tree = genome_tree)

genome_counts_filt <- genome_counts[genome_counts$genome %in% rownames(genome_gifts),]
genome_counts_filt <- genome_counts_filt %>%
  remove_rownames() %>% 
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0))%>%
  rownames_to_column(., "genome")

genome_gifts1 <- genome_gifts[rownames(genome_gifts) %in% genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]

dist <- genome_gifts1 %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, dist = dist)

7.2.1 Permanova analysis

set.seed(1234)

Richness

sample_metadata_row<- column_to_rownames(sample_metadata, "sample") 
sample_metadata_row <- sample_metadata_row[labels(beta_q0n$S), ]

betadisper(beta_q0n$S, sample_metadata_row$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df    Sum Sq    Mean Sq      F N.Perm Pr(>F)
Groups     1 0.0000047 0.00000467 0.0151    999  0.888
Residuals 22 0.0067978 0.00030899                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           Pre_grass Post_grass
Pre_grass                 0.885
Post_grass   0.90328           
adonis2(beta_q0n$S ~ diet,
        data = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))),
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
tinytable_c3mx4ut1p76rbne0e35m
term df SumOfSqs R2 statistic p.value
diet 1 0.0006002558 0.030093 0.6825872 0.785
Residual 22 0.0193464336 0.969907 NA NA
Total 23 0.0199466894 1.000000 NA NA

Neutral

betadisper(beta_q1n$S, sample_metadata_row$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.001983 0.0019829 0.1447    999  0.686
Residuals 22 0.301550 0.0137068                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           Pre_grass Post_grass
Pre_grass                 0.686
Post_grass   0.70733           
adonis2(beta_q1n$S ~ diet,
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))),
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
tinytable_i50ce727lll0g3ci24se
term df SumOfSqs R2 statistic p.value
diet 1 0.147998 0.03074149 0.6977632 0.946
Residual 22 4.666278 0.96925851 NA NA
Total 23 4.814276 1.00000000 NA NA

Phylogenetic

betadisper(beta_q1p$S, sample_metadata_row$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.001772 0.0017715 0.4816    999  0.468
Residuals 22 0.080927 0.0036785                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           Pre_grass Post_grass
Pre_grass                 0.465
Post_grass   0.49497           
adonis2(beta_q1p$S ~ diet,
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))),
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
tinytable_nyz5wunkanzv2gsc42gv
term df SumOfSqs R2 statistic p.value
diet 1 0.02764442 0.04145548 0.951464 0.59
Residual 22 0.63920164 0.95854452 NA NA
Total 23 0.66684606 1.00000000 NA NA

Functional

betadisper(beta_q1f$S, sample_metadata_row$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq     F N.Perm Pr(>F)  
Groups     1 0.008352 0.0083519 3.722    999  0.037 *
Residuals 22 0.049367 0.0022439                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           Pre_grass Post_grass
Pre_grass                  0.03
Post_grass  0.066699           
adonis2(beta_q1f$S ~ diet,
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))),
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
tinytable_t00z83wdzy7qwb6py298
term df SumOfSqs R2 statistic p.value
diet 1 0.002870065 0.03722968 0.8507252 0.452
Residual 22 0.074220722 0.96277032 NA NA
Total 23 0.077090787 1.00000000 NA NA

7.2.2 Beta diversity plots

7.2.2.1 Richness

beta_q0n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(diet) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
    geom_point(size = 4) +
    scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Post_grass" = "Captive-bred grass diet"))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(color='Origin')

7.2.2.2 Neutral

7.2.2.3 Phylogenetic

beta_q1p$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(diet) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
    geom_point(size = 4) +
    scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Post_grass" = "Captive-bred grass diet"))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(color='Origin')

7.2.2.4 Functional

beta_q1f$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(diet) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
    geom_point(size = 4) +
    scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Post_grass" = "Captive-bred grass diet"))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(color='Origin')